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CN ■ Abstract 

We present an efficient and accurate grid method for solving the time-dependent Schrodinger 
equation of atomic systems interacting with intense laser pulses. As usual, the angular part of the 
wave function is expanded in terms of spherical harmonics. Instead of the usual finite difference 






^ \ (FD) scheme, the radial coordinate is discretized using the discrete variable representation which 

^^ ' is constructed from the Coulomb wave function. For an accurate description of the ionization 

J>^' dynamics of atomic systems, the Coulomb wave function discrete variable representation (CWDVR) 
^ ; 

, P<, method needs 3-10 times less grid points than the FD method. The resultant grid points of CWDVR 



distribute unevenly so that one has finer grid near the origin and coarser one at larger distances. 

The other important advantage of the CWDVR method is that it treats the Coulomb singularity 

accurately and gives a good representation of continuum wave functions. The time propagation of 

,^£ ! the wave function is implemented using the well-known Arnoldi method. As examples, the present 



> 

oo 



o 
o 



method is applied to the multiphoton ionization of both H and H in intense laser fields. Short- 
time excitation and ionization dynamics of H by static electric fields is also investigated. For a 
<^, wide range of photon energies and laser intensities, ionization rates calculated using this method 






X 



are in excellent agreement with those from other theoretical calculations. 



This manuscript is prepared for the Journal of Chemical Physics in the Section of "Theoretical Methods 
and Algorithms " . 
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I. INTRODUCTION 

With the fast advance of modern laser technologies, lasers of various frequencies at dif- 
ferent intensities are routinely available in many laboratories [l|. Studies of the highly 
nonlinear interaction of matters with strong laser pulses have revealed many interesting 
physical and chemical processes J2| • New technologies based on these physical principles are 
under fast development and new frontiers of sciences are being opened 3[ . Newly developed 
light sources have for the first time enabled physicists and chemists to trace and image the 
electronic motion within atoms and molecules on attosecond timescale j^. 

However, multiple reactional paths and many-body nature inherent to these highly non- 
linear and ultrafast processes contribute to the complexities of theoretical interpretations 
to the experimental observations. High intensities of the applied lasers make any pertur- 
bation theories no longer applicable. It is necessary to treat the Coulomb interaction and 
the interactions from the laser fields on equal footing. Many theoretical methods have been 
designed to describe different kinds of phenomena. To mention just a few, these meth- 
ods range from strong field approximation jsf , intense- field many-body S*- matrix theory 6| , 
i?-matrix Floquet method J8[, the generalized Floquet theory 7|, time-dependent density 
functional theo ry lojl and numerical integration of the time-dependent Schrodinger equa- 



tion 



Ifl 
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, IIj, ll3| . Compared to other methods, the direct solution of the time- dependent 
Schrodinger equation (TDSE) has proved to be very versatile and fruitful in explaining and 
predicting many experimental measurements for a wide range of laser parameters. Especially 
when the laser pulse length approaches the few-cycle or sub-femtosecond regime, numerical 
solution of TDSE becomes even more appropriate and efficient. 

Nevertheless, the accurate integration of the multi-dimensional TDSE is very computa- 
tionally demanding. Even the most powerful supercomputer nowadays has only made the ab 
initio integration of the two-electron atom possible [l4|. Therefore, many approximate the- 
oretical models such as reduced-dimension models and soft Coulomb potential models have 
played an important role in understanding some of the important physical mechanisms un- 
derlying the strong field phenomena [l5[- On the other hand, it is questionable to use these 
models to simulate the realistic experiments quantitatively. For example, it was pointed out 
that the dynamical motion in three dimensions is not merely a trivial extension of what 
happens in two dimensions |ly]. It was also shown that the physics of model systems using 



the soft Coulomb potential is very sensitive to the softening parameters [17[. 

Therefore, in order to produce quantitatively correct results, it is necessary to solve TDSE 
in its full dimensions and use the real Coulomb potential with its singularity treated prop- 
erly [l8|. This is especially crucial when the singularity plays an important role in the 
problem at hand [16|, |l9|. It is the purpose of the present paper to present an accurate 
and efficient method to solve the TDSE for the multiphoton ionization of hydrogenic atom. 



Unlike the usual finite difference (FD) discretization of the radial coordinate 10|, |ll|, the 
present method discretizes the radial coordinate using the discrete variable representation 
(DVR) constructed from the positive energy Coulomb wave function. We show that the 
Coulomb wave function DVR (CWDVR) is able to treat the Coulomb singularity naturally 
and provide a good representation of continuum wave functions. 

The other advantage of CWDVR is that it needs 3-10 times fewer grid points than 
FD scheme because of the uneven distribution of the grid points, in which case one has a 
coarser grid at larger distances where Coulomb potential plays a less important role and 
wave functions oscillate less rapidly. Because the CWDVR is economic, it is a promising 
step towards a more efficient treatment of the many-electron systems which is extremely 
demanding [l^. In addition, many strong field processes, such as above threshold ionization 
(ATI), high-order harmonic generation (HHG) and dynamical stabilization etc., can be well 
understood within a single active electron (SAE) picture. Therefore, the hydrogenic atom 
serves a prototype for spherically symmetrical atomic systems interacting with intense laser 
fields. 

The rest of the paper is organized as follows. In section m we give the details of how 
we construct the CWDVR, preceded by a brief introduction to the general DVR method. 
We then present in section IIIII how we apply CWDVR to discretize the TDSE for one- 
electron atomic systems in intense laser fields. In section IIVI we present some results for 
the multiphoton ionization of H and H~ as well as the short-time dynamics of H in a static 
electric field. We show that ionization rates calculated by the present economic and accurate 
method are in excellent agreement with other theoretical calculations. Finally, we give a 
short conclusion. Atomic units (a.u.) are used unless otherwise specified. 



II. COULOMB WAVE FUNCTION DVR 

The DVR method has its origin in the transformation method devised by Harris and 
coworkers J20| for calculating matrix elements of complicated potential functions in a trun- 
cated basis set. It was further developed by Dickinson and Certain 



21| who showed the 



relationship between the transformation met 
thogonal polynomials. Light and coworkers 



lod and the Gaussian quadrature rules for or- 
first explicitly used the DVR method as a 



basis representation for quantum problems rather than just a means of evaluating Hamil- 
tonian matrix elements. Ever since then, different types of DVR methods have found wide 
applications in different fields of physical and chemical problems J23[ . There continues to be 
many efforts to construct new types of DVR and to apply DVR in the combination of other 



numerical methods 24 1 . 



Essentially, DVR is a representation whose associated basis functions are localized about 
discrete values of the coordinate. DVR greatly simplifies the evaluation of Hamiltonian 
matrix elements. The potential matrix elements are merely the evaluation of interaction 
potential at the DVR grid points and no integration is needed. The kinetic energy matrix 
elements can be calculated very simply and analytically in most cases 25|. In this section, 
we first give a short introduction to the DVR constructed from the orthogonal polynomials. 
Then we present how we construct the DVR from the Coulomb wave function, which will 
be used to solve the TDSE of atomic systems in intense laser fields in section IIIII 

A. DVR Related to Orthogonal Polynomials 



It is well known [25[ that the DVR basis functions can be constructed from any orthogonal 
polynomials, P]\[{x), defined in the domain (a, 6) with the corresponding weight function 
uj{x). Let 

VNix) = ^/uj{x)/hNPN{x), (1) 

where /iat is the normalization constant such that 

I dxVM{x)VN{x) = 5mn. (2) 

Then the cardinal function of Vn{x) is given by [2a| 

C.(-) = ^^^. (3) 



in which Xi {i = 1, 2, ..., A^) are the zeros of P/v(x), and V'j^{xi) stands for the first derivative 
oiVN{,x) at Xi. Apparently, Cj(x) satisfies the cardinahty condition 

Ci{xj) = 5,j. (4) 

One can construct the DVR basis function fi{x) from the cardinal function Ci{x) as 
follows: 

i\{x) = -^Ciix), (5) 



which, at the point Xj, gives 

fi{xj) = —=5ij. (6) 



fUJi 

We know that the integration of any function F{x) can be calculated using an appropriate 
quadrature rule associated with the zeros of the orthogonal polynomials, i.e.: 

pb N 

/ da;F(x)~^^,F(xi), (7) 

where Ui is the corresponding weight at the point Xi. From the theory of classical orthogonal 
polynomials [27[, the integration formula Eq. (jTj) is exact as long as the function F{x) can be 
expressed as a polynomial of order 2N — 1 (or lower) times the weight function uj{x). With 
the help of Eqs. (H)), Q and (0), it is easy to show that the function f*{x)fj{x) satisfies this 

condition. Therefore, the following integration can be carried out exactly: 

rb N 

/ dxf*{x)fj{x) = '^ujkf*ixk)fjixk) = 6ij. (8) 

-^^ fc=i 

As a result of Eqs. (jH)) and (jTj), any local operator V{x) has a diagonal representation in 
the DVR basis set as follows: 

rb N 

/ dxf:{x)Vix)f,ix) ^ J2''kf:{xk)V{x,)f,{xk) = V{xi)Si,. (9) 

'''^ k=l 

On the contrary, the representation of a differential operator in the DVR basis is usually full 
matrix. Nevertheless, in most cases, the matrix elements of the first and second differential 
operator 

dxf:{x)^f,{x), (10) 



dx 
and 



r d^ 

can be evaluated analytically using the quadrature rule Eqs. ((7)) or (jH)). Usually, the final 
results are very simple expressions of the zeros Xi and the number of the points A^ 



B. DVR Constructed from the Coulomb Wave Function 



Following a similar spirit, we present the DVR basis function fi(r) which is constructed 
from the Coulomb wave function. The Coulomb wave function (see, Eq. (14.1.1) of Ref. 28[) 
satisfies the differential equation 

2r/ L{L + iy 









v[r) = 0, 



(12) 



where rj is real and L is a non-negative integer. Eq. (|12p admits a regular solution F^ (r], r). 
For our purpose, we will consider the regular solution for L = and 77 < 0. Denoting 



V 



2E 



(13) 



and 



r = rV2E, 



we can write Eq. ()12|) in another form as 



dr"^ 
whose solution is given by 



v[r] 



^ 2Z' 

2E + — 
r 



v{r-) = W{r)v{r), 



(14) 



v(r) = Fo (-Z/V2E,rV2E] . 



(15) 



For any given energy E and nuclear charge Z, the solution (fT3j) has simple zeros over 
(0, 00) (see Fig. 14.3 of Ref. |28|). Similar to Eq. Q, we can define the cardinal function of 

v(r) as 



a{r) 



v{r) 



(16) 



v'{ri) r — Ti 
where r^ is the zero of v{r), and v'{ri) stands for its first derivative at rj. 

Analogy to Eq. (jSJ, one can construct the Coulomb wave function DVR basis function 



h{r) 



which, at the zero r^, equals to 



C^{r) 



U{r,) 



v{r) 



Ui v'{ri) r -Ti 



5ij. 



(17) 



(18) 



Following Schwartz J29|, Dunseath and coworkers [30| constructed an appropriate quadra- 
ture rule associated with the zeros r^ of the Coulomb wave function p5|) . By using this 
quadrature rule, the integration of a function F{r) over (0, oo) can be evaluated as follows: 

/ drF{r)^y^uJiF{ri), (19) 

where the weight Ui is given by 

oJi^^, (20) 

with 

a. = v\ri). (21) 

Using the quadrature rule Eq. (fT^. we have the orthogonality of the CWDVR basis 
functions 

/•oo ^ 

/ /;(r)/,(r)rfr^J].;,/;(r,)/,(r,)=5,,. (22) 

One can also evaluate the integration of the form (jlUp and pip using the same quadrature 
rule. We thus have 

/"°° d 

P:^ ^ J^ f:ir)-fAr)dr (23) 

N N 

^ E^'^/*(^'^)/i(^^) = E ^^^^^C7j(r,) (24) 

k=\ 

and 



^ ^it^,- 



^.. - -ll^f*(''^^2fAr)dr (25) 

AT AT 

fc=l fc=l V * J 

It is easy to show that [2^ |2l|, for the solution v{r) to Eq. (|T^ . the first and second 
derivative of the cardinal function Eq. (fTBj) are respectively given by 

C'^{n) = {l-S,,)^^—, (27) 



and 



C;(r,) = S,,^ - (1 - 5,,) ^-^— I, (28) 

dak aj (rfc - r^) 



where a^ is given by Eq. (j^U and Ck is calculated by 

Ck = akWirk), (29) 

in which W{r) = - [2E + ^] (cf. Eq. dH))). 
Finally, combining Eqs. (f^ - (PHj) gives us 

P,, = (1 - 6,,) ^^, (30) 



Tj - rj 



and 



Cj 1 

Tij = -^ijTT- + (1 - Sij) -2. (31) 

ooj [n - rj) 



with the help of Eq. (QUI) . 

Note that Eq. ()22j) is satisfied approximately for the DVR basis functions constructed 
from the Coulomb wave function, while Eq. (jH)) is satisfied exactly for the DVR basis func- 
tions constructed from the orthogonal polynomial. Nevertheless, we will show later that 
CWDVR can be applied to solve the time-dependent Schrodinger equation very efficiently 
and accurately. 

III. ONE-ELECTRON ATOMIC SYSTEM IN INTENSE LASER FIELDS 

Let us consider an effectively one-electron atomic system in a laser field. In spherical 
coordinates, the time-dependent Schrodinger equation is given by 

2^vl> (r, t) = [Hoiv) + Hj{r, t)] ^ (r, t) , (32) 

in which the field free Hamiltonian Hq is defined by 

Ho = -^V' + Vh{r) 



1 A (r-^\ - 1l^ 



Vh{r). (33) 



where V"(^(r) is the effective Coulomb potential or any kind of short (or zero) range model 
potential, which can depend on the angular momentum number /. And the orbital angular 
momentum operator U' is defined by 



sin^^^V de sin^e 



which satisfies the eigenvalue equation 

L2r,^(^,0)=/(/ + 1)1^^(^,0), (35) 

with Yim [0, (f)) being the spherical harmonics. 

In Eq. ()32j). the interaction Hamiltonian Hi{r,t) describes the interaction of the active 
electron with the applied laser pulse. In the laser parameters range of interest in the present 
work, the dipole approximation is well justified. For a linearly polarized laser along the z 
axis, Hj{T,t) is given by 

i/f ^ (r, t)=r-E{t)=r cos eE{t) , (36) 

in the length gauge, and 

Hf'\v,t) = -z-A(t) ■ V = -i-A{t)VQ, (37) 

the velocity gauge, where c is the speed of light in the vacuum. Vo denotes the 0th spherical 

I — I 

component of the gradient operator V |31;j. The electric field strength E(t) is related to the 
vector potential A(t) of the laser pulse by 

E(t) = -^|A(t). (38) 

A. Discretization of the Spatial Coordinates 



In order solve the TDSE (jHzj) . we need to discretize this equation. We expand the angular 
part of the wave function in terms of spherical harmonics. The radial coordinate r can be 
discretized in different ways. The most straightforward one is to use the finite difference 
method, in which case the first and second derivative with respect to r in Eq. ()33|1 are 
approximated by formulas involving only several neighboring points. In the present work, 
we expand the radial part of the wave function in CWDVR basis functions constructed 
above. As seen from Eqs.( IHUj) and ( IHT|) . CWDVR is a global method, in which case the 
representations of the derivatives involve all the grid points. 

1. Expansion of the Angular Part 

Expanding the time dependent wave function in the following form 

'flrnir, t) . 

r 



^ (r, t)^^> (r, eA.t) = Y.Y. ^^^^^^^Yim {e, <P) , (39) 



/=0 m=-L 



and substituting into Eq. (|52|). we arrive at 

L L 






/=0 m=-L 

t t >-(» 

Z=0 m=-L 

L L ^ 

+ Y.11 Hj{r,t)-^im{r,t)Y, 






^im{r,t) 



(40) 



i=0 m=-L 

with the help of Eq. IH^ 

Multiplying the above equation by Y^^^, [6, (p) jr and integrating over r^ sin OdOdcf) on both 
sides, we can rewrite it as 



dt 



where we have made use of 



d 1 (f 



2dr2 



(41) 



2n 



sin 9d9 / d4)Yi*^, {9, 0) Yim {9, 0) = SvAn'm- 



In Eq. fl41|l . we have defined an effective potential term 



(42) 



VefAr) ^ V^{r 



l{l + l) 
2r2 



(43) 



and a laser-interaction-related term, 



L L 



^ -^ poo PIT /•27r 

«=0 m=-L "^0 -^0 -^0 

X ^Y;^, {9, 0) Hj{t, t) ^-^Ur, t)Yi^ {9, 0) . 

In the length gauge, substituting Eq. PUJ) into Eq. (jUj) and making use of Eq. 
the following formula [^ 



cos 9Yi 



Im 



O'l+lmYl+lm {9, 0) + aimYl-i 



m v^5 r; ) 



where 



air. 



(I — m) {I + m) 
{21 -I) {21 + 1)' 



(44) 
and 

(45) 
(46) 
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we arrive at 



r(i) 



Hr{r,t) 



U'm' 



rE{t) [avm'^i'-im'{r,t) + av+im'^v+im'{r,t)] . 



In the velocity gauge, substituting Eq. (jTfj) into Eq. (jUJ and making use of Eq. 
the formula 



-R{r)Yi 



Im 



ai+imYi+im. {6,(j))- { — — ] R{r) 



we arrive at 



r(^) 



Hr\r,t) 



U'm' 



r \dr r 

+aimYi-im (^, 0) - ( :r + - 1 ^(^) 
r \ar r 



%A{t)- [rai'jn'Vl'^lm'{r,t) - (/' + l)aif+im'^l'+lm'{r,t)] 



-iA{t)— [ai>m'^i'~im'{r,t) + ai'+im'(pi'+im'{r,t)] 



with aim given by Eq. pH|l . 



2. Discretization of the Radial Coordinate 



(47) 
and 



(48) 



(49) 



We have now changed the TDSE (J22I) into ()41|) with the interaction term [Hi{r,t)]f,^, 
given by Eq. ()47|) in the length gauge and by Eq. ()49p in the velocity gauge. Although we 
deal with the r coordinate using CWDVR in the present work, we still give the FD formulas 
for the purpose of comparison. In the FD scheme, the grid points are chosen to be equal 
spacing Ar such as 

ri = iAr, i = l,2,...,N. (50) 

The 5-point central finite difference approximation to the first and second derivative of 
<fim{r,t) are given by 



—(pim{r, t) = — — — [ipim{r - 2Ar, t) - 8ipim{r - Ar, t) 
ar IzAr 

+8ifim{r + Ar, t) - ipimir + 2Ar, t)] , 



(51) 



and 



ar"^ 



[^im{r - 2Ar, t) - 16(/)/m(r - Ar, t) 



12{Arf 
+30<f im{r, t) - IQ^imir + Ar, t) + i^im{r + 2Ar, t)] 



(52) 
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In order to have give a better ground state energy, we have used the following approxima- 
tion for the second derivative at the first point to take into account the boundary condition: 

d^ 1 

dr^ 12 (Ar) 

+Coipimir,t), (53) 



-J3}Pim{r,t) = -^^:^y^^-^[30ipim{r,t) -16ipim{r + Ar,t) + Lpim{r + 2Ar,t)] 



where Cq is a constant which depends on the grid spacing Ar. For example, we take 
Co = —1.48986 for Ar = 0.2 and get a converged H ground state energy of —0.500000065 
a.u., but Co = —1.814116 for Ar = 0.3 in order to get the same value of the ground state 
energy. 

Apparently, one needs a very small spacing Ar in order to have a good approximation 
to these derivatives. For the atomic systems interacting with the intense laser pulses, the 
electronic wave function can be generally driven to hundreds of or even thousands of atomic 
unit away from the nuclear. Therefore, we need a very large number of the grid points. 
Another disadvantage of the FD scheme is that, one needs to carefully deal with the Coulomb 



singularity at the origin 



0. 

Now let us turn to discretize r using the CWDVR basis functions which we discussed in 
the last section. As we mentioned before, the CWDVR has several advantages: firstly, it 
deals with the singularity of the Coulomb-type potential at the r = naturally; secondly, 
the grid points {i.e., zeros of Coulomb wave function) are dense nearly the origin where the 
Coulomb potential plays a crucial role and sparse at large distances where it is not very 
important; thirdly, compared to FD scheme, much fewer grid points are needed for the same 
extent of the grid because of the uneven distribution of the grid points. 

The CWDVR basis functions /^(r) are given in Eqs. dlTD , dUD and (J22D. Let us start 
from Eq. (jlT|l and expand ipim{r,t) in terms of /i(r) as follows: 

N 

ipi^{r, t) = Y^ Dum{t)fi{r), (54) 

i=l 

where the coefficient is given by, 

Dum{t) = / drf*{r)(pim{r,t) 
Jo 

N 

fc=i 
= yMiV^im{ri,t). (55) 

12 



Substituting Eq. (|^ into Eq. (HH), multiplying both sides by f*,{r) and integrating over 
r, we arrive at 

N 

dt 

+ \HAr,t)]^n'^, (56) 



d 

i-^Difi'm'it) = 2_^ TifiDiitm'it) + Veff (rj/) Di'i'm'it) 
i=l 



where we have made use of Eq. (J2SI)- [Hi{t)]'^,i,^, stands for the matrix element of the 
interaction term, which is shown to be 

Hi^\t) =52/ drfif{r)fi{r)E{t)r[ai-,n'Dii'-im'{t) + ai'+im'Dii>+i^>{t)] 

L J i'l'm' ^ Jq 

= E{t)ri, [ai'm'Dii'.im'it) + «//+!„/ AF+lm'(t)] • (57) 

in the length gauge. In the velocity gauge, it takes a slightly complicated form as 

-lip ^ r°° 1 

Hi^\t) = iA{t) V / drfi\r)-U{r) [l'av^,Du'^im'{t) - (/' + 1) ar+wAF+w(t)] 

-iA{t) ^ / drfi,{r)—fi{r) [ai'm'Du'^im'{t) + ai'+im'A«'+im'(^)] 



i=i-^o 



iA{t) [l'aiirn'Diq/_im'{t) — (/' + 1) Oi'+lm' A'i'+lm' (^)] 



N 



—iA{t) 2_^ Pi'i [(i'l'm' Dill -im' if) + O'l'+lm'Diii+im'it)] . (58) 

i=l 

where we have made use of Eq. ((221) • The matrix elements of Pin in Eq. (J3H|) and Tin in 
Eq. (jSni) are calculated analytically from Eqs. (jHSI) and (p?T|l respectively. The zeros r^ needed 
for evaluating these matrix elements are calculated with the help of COULFG J22l- Please 
note that, for a linearly polarized laser considered in the present work, the subscript index 
m' is equal to 0. In this case, we only have two dimensional matrices with indexes i' and /'. 

3. Distribution of the CWDVR Grid Points 

The grid point r^ of CWDVR is the solution of v{r) = where v{r) is defined by Eq. (|15p 
for any given energy E and nuclear charge Z. Therefore, the distribution of the CWDVR 
grid points can be adjusted by the value of parameter Z and k, = \plE. 

In tablelU we compare the grid points distribution of the CWDVR for the same maximum 
grid point value r^ax — 150 and for different k values ranging from 0.5 to 5. The parameter 
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TABLE I: Comparison of the 


grid points 


distribution of CWDVR for different Z and k. 




K = 0.5 


K = 1.0 


K = 2.0 


K = 3.0 


/€ = 4.0 


«; = 5.0 


Z=12 


N 


47 


64 


106 


152 


198 


245 




n 


0.1529 


0.1526 


0.1517 


0.1501 


0.1480 


0.1455 




TN 


153.29 


150.70 


150.04 


150.86 


150.45 


150.52 




Aro 


0.3589 


0.3565 


0.3472 


0.3335 


0.3171 


0.2994 




Ari5 


2.3650 


1.9090 


1.3227 


0.9655 


0.7491 


0.6091 




Ari5o 


4.9116 


2.9158 


1.5402 


1.0381 


0.7815 


0.6263 


Z=20 


N 


56 


72 


112 


156 


202 


248 




n 


0.09174 


0.09169 


0.09148 


0.09114 


0.09066 


0.09007 




TN 


150.27 


151.39 


150.74 


150.47 


150.73 


150.43 




Aro 


0.2157 


0.2151 


0.2130 


0.2097 


0.2053 


0.2001 




Ari5 


1.8308 


1.6554 


1.2137 


0.9132 


0.7270 


0.5972 




Ari5o 


4.3561 


2.7914 


1.5209 


1.0320 


0.7789 


0.6250 



Z takes either 12 or 20. We list in this table the number of the grid points A^ (up to the 
first point which is greater than 150) and the value of the first and the last grid point ri and 
ttv- We also give three grid spacings Aro, Aris and Ariso which correspond to the spacing 
between the first two points, the two points around 15 and the last two grid points. 

One notices that the value of Z mainly determines the value of the first grid point ri, 
i.e., the greater Z is, the smaller ri becomes. However, increasing the value of k will mainly 
decrease the spacing between the grid points at large r and thus one needs more grid points 
for the same value of r^ax = 150. At the same time, larger k means much more even 
distribution at larger distances. For example, Aris and Ariso for k = 5 differ much less 
than that for k, = 0.5. 

For the same Vmax = 150 a.u., one observes that the number of the grid points for CWDVR 
cases is only 1/10 to 1/3 of that for FD case if we choose even spacing Ar = 0.2 a.u.. As 
we will discuss later, most of our results for the hydrogen are convergent for k, = 1.0 and 
Z = 20 case where only 1/10 of the number of the grid points for FD is needed (A^ = 72 vs 
N = 750). 
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B. Wave Function Propagation in Time 

After the discretization of the spatial coordinates, one has to advance the evolution of 
the initial wave function in discretized time. If the initial wave function one desires is the 
electronic ground state, it can be computed by the propagation of any trial wave function in 
imaginary time. The Schrodinger equation in imaginary time is mapped into the diffusion 
equation. In this case, any excited states will decay away faster than the ground state 
since the latter has the lowest energy. Once the energy is adjusted to the true ground 
state energy Eq, the asymptotic solution is a steady-state solution. One will thus obtain a 
converged ground state energy and wave function on the spatial grids after a sufficiently 
long time of diffusion. 

Concerning the propagation of the time-dependent Schrodinger equation, there exist 



many differen 
Taylor series 



Q, 



Chebychev polynomial expansion 
There have been a number of authors 



] methods such as split operator 
and Arnoldi/Lanczos method 
who made a detailed comparison study of the efficiency and accuracy among different propa- 
gation schemes J37|. Although the choice of the propagation scheme depends on the internal 
characteristics of the physical problem at hand, it is generally accepted that Arnoldi/Lanczos 
method proves to offer an accurate and flexible approximation of the matrix exponentials 
involved in the propagation of wave function for most practical applications (Note that 
Arnoldi method can apply to nonsymmetric matrices, but Lanczos method only applies to 
symmetric or Hermitian matrices). 

In the Arnoldi/Lanczos method, the eigenvalue and eigenstate of a large matrix is ap- 
proximated in the Krylov subspace spanned by the vectors {tpo, HtpQ, H'^tpQ, ..., H^^ipo} for 
the order m which is typically much smaller than the dimension of the Hamiltonian ma- 
trix H itself. In practice, we need to use Arnoldi/Lanczos algorithm J38[ to generate m 
orthonormal basis vectors Kn+i in the Krylov subspace. At the same time, an upper Hes- 
senberg matrix /im+i is formed resulting from these processes. The eigenvalues of the original 
large matrix are thus approximated by the eigenvalues of a Hessenberg matrix of dimension 
772-1-1. Therefore, Arnoldi/Lanczos method is especially suitable for the large matrix and 
multi-dimensional problems. 

An extensive software package "Expokit" for the computation of matrix exponential was 
developed by Sidje j^. In this package, the author has provided several alternatives to 
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compute the matrix-valued function e*^ for small, dense complex matrices H. In addition, 
"Expokit" has functions for computing e^^tpo for both small, dense matrices and large, sparse 
matrices. 

In the present work, we have incorporated this software into our code for the accurate 
propagation of the wave function for atomic systems interacting with strong laser pulses. In 
each time step, we first generate the orthonormal vectors Vm+i and the upper Hessenberg 
matrix hm+i for the wave function \E'(t) using the Arnoldi process. Then the exponential 
of the upper Hessenberg matrix e~*'^*'^™+i is carried out by the irreducible rational Pade 
approximation combined with scaling-and-squaring. Finally, the wave function at the next 
time step is given by 

^{t + At) = \^{t)\V^+ie-'^"'^+'ei, (59) 

where |\l^(t)| stands for the norm of the wave function at time t and ei is the first unit basis 
vector J39[. 

For all the results presented in this paper, we use the Arnoldi order m = 30 and the 
propagation time step At = 0.01 a.u.. However, we find that for some of the laser parameters 
we consider, the results are already converged for much lower order of m (e.g., 12) and much 
larger time step At (e.g., 0.04 a.u.). 

IV. RESULTS AND DISCUSSIONS 

In this section, we will present some numerical results for ionization of the hydrogen 
atom by strong static electric fields and by intense laser fields. In order to illustrate that 
our method can equally apply to any effectively one-electron atomic system with a proper 
SAE model potential, we have also calculated multiphoton detachment rates of the negative 
hydrogen ion for several laser frequencies at different intensities. Excellent agreements are 
achieved between our results and the previous theoretical calculations. 

For an accurate evaluation of the ionization rate for a particular laser intensity, the laser 
pulse should ramp on adiabatically and keep constant for a sufficiently long time. The laser 
pulse ramping-on time must be large compared to the initial atomic orbital period of the 
bound electron and thus the electron energy will adjust adiabatically to the rising intensity. 
The constant-intensity time should also be sufficiently long so that the frequency bandwidth 
is small compared with the laser frequency. Only in this case, we can treat the laser field as 
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monochromatic. Taking these into consideration, we take the vector potential A(t) of the 
laser pulse having the following form: 

A(t) = Aof{t) cos{ujt)k, (60) 

where the pulse envelope, f{t), is given by 



1 



1 — cos 



nt 



f{t) = { n v-yj' - - ^' (61) 

1, n < t < Ti + T2, 

where ri and T2 are normally taken to be 5-10 and 10-40 laser cycles, depending on the laser 
frequency under consideration. The peak value of the vector potential Aq is connected to 
the peak laser intensity Iq by 

Ao = ^ = -^/f , (62) 

UJ UJ\I lau 

where the atomic unit of laser intensity lau equals to 3.5094 x 10^^ W/cm^. The electric 
field strength E(t) is calculated from Eq. (jHHj) . 

In order to avoid the reflection of wave function by the edge of the r grid, the wave 
function is multiplied by an absorbing function in each step of the propagation as follows 

^(r,t) =M(r)^(r,t), (63) 

in which 




r < r, 



at 



(64) 



where r^ = armax and To- = aVmax with r^ax being the maximum value of the grid. It is 
very important to carefully choose the absorbing parameters a and a such that the function 
M{r) is sufficiently smooth and that the wave function near the edge is completely absorbed 
without any reflection. It is extremely crucial to avoid any unphysical effects induced by 
the absorbing potential J40|. According to our experiences, we take 0.3 < a < 0.6 and 
0.5 < cr < 5.0 depending on the value of r^ax and the laser parameters. The criteria is that 
the time-dependent physical quantities, such as the population decay of the ground state, 
should be converged against any small variation of these parameters in the vicinities of some 
values. 

Of course, convergence has to be achieved against changing of other parameters such as 
the largest quantum number of the angular momentum L, the number of grid points N in 
r coordinate, the propagation time step At and the Arnoldi order m etc.. 

17 



0.0 



-0.5 



-1.0 



-1.5 



g 

Q. 
O 
CL 



O) -1.2 
O 



_ FD, Ar=0.2 

K=0.5, Z=20 




(a) - 


- K=1.0, Z=12 






K=1.0, Z=20 




^*=^te~^ 


K=2.0,Z=12 




^'^^^^^N 


K=2.0, Z=20 
1 


1 1 1 


1 , , , , ''" 



10 



15 



-1.4 



_ 


-_ 1 

— ^ ^ 




1 




1 


(b) 


_ 


- 




- 





■--- 


■ V . 





- 








' '"»>-,. ""-* 


- 


1 




1 




1 




_ 







14.1 



14.4 14.7 

Time (Cycle) 



15.0 



FIG. 1: (color online). Convergence of the CWDVR grid for different k and Z values. The natural 
logarithm of the population decay as a function of time (in unit of laser cycle) are shown for 
the laser wavelength 780 nm at the peak intensity of 2 x lO"*^^ W/cm^. (a) Results from different 
CWDVR grids {k and Z as indicated) are compared against that from FD method; (b) A magnified 
version of (a) for the time from 14 to 15 cycle. 



A. Choice of Gauge and Convergence of the CWDVR Grid 

Although different gauges describing the interacting Hamiltonian Hj{t) should in principle 
lead to the same physical results, it is not true in practice because of different approximations 
and inaccurate numerical wave functions. For some approximate methods, length gauge 
proves to be better than the velocity gauge in some circumstances 41| . However, as discussed 
by Cormier and Lambropoulos 42| , velocity gauge is preferable for some ab initio numerical 
calculations. This is connected with the fact that the canonical momentum in the velocity 
gauge has been reduced to a slowly varying variable since the momentum due to the strong 
field has been mainly removed 42|. In this case, we can avoid some widely varying variables. 
Especially in the present case where the wave function is expanded in terms of the spherical 
harmonics, much smaller L is needed for a converged result in the velocity gauge than that 
needed in the length gauge. 

In the present work, we have performed our results in the velocity gauge for H in laser 
fields. However, length gauge is used for H in static field and H~ in laser fields. One of 
the reasons for the latter choice is that it is easy for us to compare our results with the 
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FIG. 2: The logarithm of the decay of the ground state and population within different spheres of 
r for H at the laser frequency lo = 0.6 a.u. and the peak intensity /q = 4.375 x 10^'^ W/cm^. 

previous results. Another reason, for H~, is that the angular- momentum-dependent model 
potential |43| makes the form of the potential in the velocity gauge is not able to be obtained 
obviously J4^. 

Since the CWDVR grid depends on the two parameters Z and k, one has to make sure 
that any calculation should be converged against changes of them. As shown in table U 
the CWDVR grid for smaller k is much coarser. At the same time, the first point of the 
grid is mainly decided by the value of Z. It is expected that much finer CWDVR grid will 
have a better representation of the Coulomb potential and the wildly driven electronic wave 
function by the laser field. 

In order to test the convergence of CWDVR grid, the multiphoton ionization of H for 
an infrared laser at a moderate laser intensity serves as a good example. For a hydrogenic 
atom, the potential is given by 



Vc(r) = 



(65) 



where Z is the nuclear charge and equals to 1 for H. For the CWDVR grids we consider in 
table m we get the H ground state energy of —0.49999997 a.u. or better. 

In figure n we present the natural logarithm of the population decay within a sphere of 25 
a.u. of r calculated by the FD method and by the CWDVR grids corresponding to different 
Z and K. With the wavelength A = 780 nm and peak intensity Jq = 2 x 10^^ W/cm^, the 
laser pusle has a 5-cycle ramping-on time and keeps constant for another 10 cycles. The 
maximum number of angular momentum L is taken to be 20 for a converged calculation. In 
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TABLE II: Ionization rate T (in s ^) for ionization of H by a linearly polarized laser of intensity 
Iq (in W/cm^) and frequency lo fin a.u.) . The present results (a) are compared with the previous 



results: (b) Chu and Cooper 



43]; (c) Pont, Proulx and Shakeshaft J47|; (d) Kulander |10|. The 



Id; 



number with parenthesis like p(q) is understood as p x 10'^. 



UJ 


/o 


pa pb pc pd 


0.55 


7.00(12) 


1.43(13) 1.43(13) 1.43(13) 1.4(13) 


0.28 


7.00(12) 
4.38(13) 


3.73(11) 3.73(11) 4.0(11) 3.3(11) 
1.33(13) 1.33(13) 1.35(13) 1.2(13) 


0.20 


4.38(13) 
1.75(14) 
3.94(14) 


3.74(12) 3.86(12) 4.0(12) 2.8(12) 
2.65(14) 2.89(14) 2.7(14) 4.0(14) 
6.14(14) 5.65(14) 6.0(14) 7.0(14) 



the FD case, the grid spacing Ar is taken to be 0.2 a.u., and the number of the grid point 
A^ = 750 for the maximum value of the grid point r^ax = 150. The corresponding number of 
the r grid points used for CWDVR cases are indicated in table HI From figure ^ we conclude 
that our results are indeed converged to the FD difference results as k is increased from 0.5 
to 2.0. Even for the coarsest case for k, = 0.5, in which case only 56 grid points are used, 
the result is reasonably good. The ionization rate estimated from this curve is 4.63x10^^ 
s~^ which is close to the FD result 5.07x10^^ s~^ and the CWDVR result 5.05x10^^ s~^ for 
K = 2. We also observe that, for the present case, results are not very sensitive to changing 
of the value of Z{ioi the same value of n). However, we find that, in most of the cases 
that we consider later, Z = 20 is usually more preferable for its better representation of the 
Coulomb potential near r = 0. 

In practice, the ionization rate is fitted by an exponential decay of the ground state 
population or the decay of the total population within some sphere of r. In the present 
calculations, we have defined an inner sphere and a middle sphere with radius rj„„ = 25 and 
fmid = 50 a.u. respectively. The population remaining within these two spheres Pj„„ and 
Pmid, together with the probability in the ground state Pgrd and the total population in the 
whole box Pwhh are recorded for each time time step. As an example, we show in figure |21 
for the ionization of H at the laser frequency 0.6 a.u and peak intensity 4.375 x 10^^ W/cm^. 
The laser has a ramping on time of 5 cycles and keeps constant for 50 cycles. We observe 
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TABLE III: Multiphoton ionization rates of H at different laser electric field strengths Frms for 
different photon energies uj. Results calculated from the present method are compared with those 
of Chu and Cooper by a nonperturbative L^ non-Hermitian Floquet method 45 1. 











r/2 


(a.u.) 








u> 


p — 

^ rm.s — 


0.01 a.u. 


Frms = 0.025 a.u. 


7? — 

^ rm.s — 


0.05 a.u. 


Frms = 0.075 a.u. 


(a.u.) 


Present 


Ref [45] 


Present 


Ref [45] 


Present 


Ref [45] 


Present 


Ref [45] 


0.60 


0.125(-3) 


0.125(-3) 


0.783(-3) 


0.784(-3) 


0.313(-2) 


0.314(-2) 


0.704(-2) 


0.711(-2) 


0.55 


0.173(-3) 


0.173(-3) 


0.108(-2) 


0.108(-2) 


0.435(-2) 


0.436(-2) 


0.982(-2) 


0.989(-2) 


0.50 


0.250(-3) 


0.247(-3) 


0.157(-2) 


0.154(-2) 


0.647(-2) 


0.624(-2) 


0.149(-1) 


0.139(-1) 


0.30 


0.378(-5) 


0.377(-5) 


0.131(-3) 


0.131(-3) 


0.160(-2) 


0.161(-2) 


0.594(-2) 


0.639(-2) 


0.28 


0.451(-5) 


0.451(-5) 


0.161(-3) 


0.161(-3) 


0.204(-2) 


0.204(-2) 


0.821(-2) 


0.815(-2) 


0.27 


0.503(-5) 


0.502(-5) 


0.180(-3) 


0.180(-3) 


0.230(-2) 


0.231(-2) 


0.920(-2) 


0.920(-2) 


0.26 


0.562(-5) 


0.562(-5) 


0.202(-3) 


0.202(-3) 


0.256(-2) 


0.261(-2) 


0.106(-1) 


O.llO(-l) 


0.22 


0.335(-6) 


0.180(-6) 


0.200(-4) 


0.189(-4) 


0.474(-3) 


0.511(-3) 


0.109(-1) 


0.173(-1) 


0.20 


0.140(-6) 


0.106(-6) 


0.452(-4) 


0.467(-4) 


0.321(-2) 


0.350(-2) 


0.743(-2) 


0.683(-2) 


0.19 


0.180(-4) 


0.715(-5) 


0.454(-3) 


0.446(-3) 


0.213(-2) 


0.214(-2) 


0.466(-2) 


0.437(-2) 


0.18 


0.884(-6) 


0.791(-6) 


0.945(-4) 


0.948(-4) 


0.158(-2) 


0.148(-2) 


0.243(-2) 


0.376(-2) 



that the population decay within different spheres of r are parallel lines with the line of the 
ground state decay. Therefore, it does not matter in this case which curve to be used for the 
estimation of the ionization rate. The calculated ionization 1.5658x10"^ a.u. in this way is 
in very good agreement with Chu and Cooper's result 1.5672x10^'^ a.u. (45| . 

In a similar way, we have also estimated the multiphoton ionization rate of H by the 
laser of wavelength 1064 nm at the peak intensity 1 x 10^^ W/cm^. The result is 2.97x10^^ 
s~^, which is in good agreement with an independent molecular code result 2.85x10^^ s~^ 

n 

in cylindrical coordinates and the FD result 2.92x10^^ s ^ in spherical coordinates J46|. 

B. Multiphoton Ionization of H by Intense Laser Pulses 

From the above test calculations, we have seen that the present method is applicable to 
different laser wave lengths at different intensities. In this section, we give more rigorous 
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FIG. 3: (color online). Depletion of the ground state of H by a static electric field with the field 
strength 0.08 a.u.. The ground state probability is shown as a function of the field duration. 
Results calculated by different CDVR grids for different k values are compared against the result 
calculated by FD method with Ar = 0.1 a.u.. 

tests of the present method by calculating the ionization rates of H for a large number of 
wavelengths from very low to very high peak intensities. 

For all the calculations presented here, we use the CWDVR grid for k = 1.0 and Z = 20. 
Therefore the number of grid point A^ is 72 with r^ax = 151.39 a.u.. The maximum number 
of angular momentum L = 10 gives converged results in all the cases for the velocity gauge 
description of the interaction term. The parameters for the absorbing function (jMj) a and 
a are taken to be 0.4 and 4 respectively. In table m we compare ionization rates calculated 
using the present method with other theoretical results. We notice that our results agree 
better with those of Chu and Cooper J4^ by an ab initio Floquet calculations. However, the 
present time-dependent method results are also in reasonable agreement with those TDSE 
calculations of Ref. |47[ using a complex Sturmian basis and those of Ref. jlO| using a FD 
method. The results oi u = 0.2 a.u. for intensities 1.75 x 10^^ and 3.94 x 10^^ W/cm^ are of 
less favorable agreement because high ionization rates lead to very fast decay of the ground 
population and thus estimations from TDSE methods become less accurate. 



In table 
Cooper 



TTTj we compare ours results with the benchmark calculations done by Chu and 
using the ab initio nonperturbative L^ non-Hermitian Floquet method. The 



laser frequency u vary from 0.6 to 0.18, which correspond to one- to three-photon ionization 
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FIG. 4: H ground state survival probability Pis as a function of time in static electric fields. The 
field strength F is taken to be: (a) 0.005; (b) 0.04; (c) 0.06; (d) 0.08 a.u.. These results are in 
perfect agreement with those by Durand and Paidarova |5l| and by Scrinzi [49|. 



of the ground state H. Please note that Eq is connected with Frms by 



Eq = V2Frms = yl/h 



Again, very good agreement are achieved for these wide ranges of laser parameters. 



(66) 



C. Ionization of H by Static Electric Fields 

Now let us turn to the ionization of H from the ground state by a static electric field. 
Excitation and ionization dynamics of atoms and molecules by static electric fields has been 



of great importance since the foundation of quantum mechanics 



great interest of current study 



49, 



50 



4S I and continues to be of 



Although the long time behavior is dominated by the exponential decay of the ground 
state of H, large deviations from the naive exponential decay are expected because of the 
sudden turn-on of static electric field. Actually, the spectrum of the Hamiltonian for this 
system is unbounded from below. There are many resonances present. The deviation from 
exponential decay is expected to be large for strong fields cases JSIJ. However, there is a 
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substantially long time of transitory region for weak and intermediate fields as well, whose 
length depends on the particular field strength. 

In order to investigate this transitory regime, one has to make sure that the nonexpo- 
nential decay indeed comes from physical dynamics rather than numerical reasons such as 
nonconvergence of the grid or reflection of the wave function from the edge. This point is 
extremely important for the static electric field case where the wave packet is driven away 
in one direction rather than in a oscillatory fashion as for a laser pulse case. In figure El we 
show the logarithm of the ground state probability of H by a static electric field F = 0.08 
a.u. for different CWDVR grids. Note that we use the length gauge in this case for the 
dipole interaction. The maximum of the grid r^ax for all the cases is taken to be around 150 
a.u. and the corresponding number of the grid points are listed in table Q] We take Z = 20 
for all the CWDVR grids. In addition, the absorbing function parameters a and a are taken 
to very extreme values of 0.25 and 0.4 respectively in order to avoid the reflection from the 
edge. For this field strength, we can get converged result only if k > 4.0, which case the grid 
is dense enough for a good representation of the interaction term using the length gauge. 
Note that the performance of the present CWDVR is much better than the DVR method 

n 

used by Dimitrovski and coworkers [52| , in which case the converged result is only achieved 
for t ^ 40. Therefore, one has be be very careful not to interpret the nonexponential decay 
for K < 3.0 as any physical excitation and ionization dynamics J5l|. 

Using the CWDVR k > 4.0 and Z = 20, we have studied the short-time dynamics at 
other field strengths. In figure IH we present the population decay of the ground state for 
different electric field strength F. For the lowest field strength F = 0.005 a.u., one observes 
both a small- and large-time scale regular oscillation, which implies reversible processes. 
For the case when F becomes 8 times of that in (a), we observe a quadratic decay in the 
first 8 a.u., followed by a irregular transition time before the system becomes stable after 
90 a.u. or so (it will be a very slow exponential decay afterwards, with F ~ 10~^ a.u.). As 
the field strength increases further in (c) and (d), the transitory time becomes shorter and 
shorter, which is followed by a purely exponential decay. The transitory regions are shown 
by Durand and Paidarova [5l| to directly related to the 2s-2p resonances by inspecting the 
spectral density line-shape. Note that, the fast quadratic decay is present for (c) and (d) as 
well when the field is turned on. 

It is remarkable to note that the present time-dependent calculations for the entire region 
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TABLE IV: Ionization rate T (in a.u.) for ionization of the ground state of H by a static electric 



field of strength F (in a.u.). Results are compared with: (a) Scrinzi 



Bauer and Musler 



53|. 



4g|; (b) Peng et al [i$; (c) 



Li; 



F 



Present 
Ref [42] 
Others 



0.06 



0.08 



0.1 



0.5 



5.1509(-4) 4.5396(-3) 1.42(-2) 5.64(-l) 

5.1508(-4) 4.5397(-3) 1.45(-2) 5.60(-l) 

5.15(-4)'' 4.55(-3)^ 1.2(-2)^ 5.4(-l)^ 



from very weak to very strong field strengths are in complete agreement with the results 
using the complex scaling methods (49l . l5ll |. This is further confirmed by comparing our 
ionization rates fitted by an exponential decay in the region where the transitory time is 
over, which is shown in tabl e IIVI Also shown in this table are the results of some other 
time-dependent calculations |l3|, |53| . 



D. Multiphoton detachment rates of H by strong laser pulse 

In order to show that the present method can equally apply to any atomic systems with 
appropriate SAE model potential, we study the in this section the multiphoton detachment 
of the negative ion of hydrogen. For H~,_we use the angular- momentum-dependent model 
potential proposed by Laughlin and Chu 



1 



2r Oid 



Vh{r) = 1 + - e"^'- - -^W, {r/r,) + n,(r). 



2r4 



where 



Wj{x) 



(67) 



(68) 



ui{r) = (co + cir + C2r^) e'^^ (69) 

In the present calculations, we use the same values of parameters listed in Table I. of Ref. J43| . 
The length gauge is used to describe the interaction term. Please also note that, we adopt 
the corresponding dipole operator which includes the contribution induced on the hydrogen- 
atom core by the outer loosely bounded electron. Instead of r, the dipole operator D in this 
case is given by 



D 



25 



Ws {r/rc) 



(70) 



TABLE V: Multiphoton detachment rates of H for 1064 nni, 1640 nni and 1908 nni at different 
intensities ranging from 1 x lO^*' to 1 x 10^^ W/cm^. The present results from the CWDVR are 



compared with those of Haritos and coworkers 



5g] and those of Telnov and Chu 



57 



Sg]. The 



detachment rate is quoted as p{q) which stands for p x 10'^. 





Photodetachment Rate (a.u.) 


Intensity 


Present 


1064 nm 
Ref [^ 


Ref [^ 




1640 nm 






1908 nm 




(W/cm2) 


Present 


Ref [^ 


Ref [^ 


Preseni 


Ref [^ 


Ref [^ 


1 X 10^° 


4.56(-5) 


4.50(-5) 


4.65(-5) 


6.5(-7) 


5.05(-7) 


2.98(-7) 


4.90(-7; 


4.33(-7) 


4.80(-7) 


5 X 10^° 


2.25(-4) 


2.20(-4) 




7.7(-6) 


6.82(-6) 




1.14(-5; 


1.04(-5) 




8 X 10^° 


3.57(-4) 


3.56(-4) 




1.81(-5) 


1.71(-5) 




2.81(-5; 


2.58(-5) 




1 X 10^^ 


4.43(-4) 


4.43(-4) 


4.53(-4) 


2.79(-5) 


2.63(-5) 


2.78(-5) 


4.27(-5; 


3.98(-5) 


4.37(-5) 


2 X 10^1 


8.58(-4) 


8. 70 (-4) 




1.01(-4) 


0.99(-4) 


1.05(-4) 


1.55(-4; 


1.46(-4) 


1.58(-4) 


3 X 10^1 


1.25(-3) 


1.28(-3) 




2.14(-4) 


2.10(-4) 




3.23(-4; 


3.01(-4) 




4 X IQii 


1.61(-3) 


1.66(-3) 




3.68(-4) 


3.55(-4) 


3.76(-4) 


5.18(-4; 


4.95(-4) 


5.30(-4) 


6 X 10^^ 


2.24(-3) 


2.38(-3) 




7.62(-4) 


7.19(-4) 




9.52(-4; 


9.78(-4) 




8 X 10^^ 


2.66(-3) 


2.98(-3) 




1.16(-3) 


1.15(-3) 


1.23(-3) 


1.48(-3; 


1.55(-3) 


1.52(-3) 


9 X 10^^ 


2.81(-3) 


3.24(-3) 




1.39(-3) 


1.40(-3) 




1.72(-3; 


1.88(-3) 




1 X 10^2 


3.05(-3) 


3.46(-3) 


2.95(-3) 


1.64(-3) 


1.64(-3) 


1.73(-3) 


2.16(-3; 


2.12(-3) 


2.18(-3) 


Use th 


is model 


potential, 


the ground state 


we get is 


-0.027730 a.u. for all the CWDVR 


grids liste 


d in tabl 


elU which 


is in very good ag 


reement with the value —0.027733 a.u. calcu- 



lated by Telnov and Chu pj| and with the experimental measurement —0.027716 a.u. |55l |. 

In table I3 we compare the total detachment rates of H~ calculated from the present 
method with other theoretical calculations. We get converged results of the CWDVR grid 
for K, = 2 and Z = 20. It is not surprising that the requirement of the grid is less demanding 
for H^ case than for the H in static electric field case since we have a short-range potential 
here. For all the laser parameters considered in table I3 we get very good agreement with 
those results of Refs. j^ and j^, |^ . 

However, Telnov and Chu use the dipole moment r, which accounts for the slight differ- 
ences between their results and ours. We have tested our code with the latter dipole moment. 
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we get detachment rate to be 4.51 x 10"'* a.u. instead of 4.43 x 10"'* a.u. for 1064 nm at 
1 X 10^^ W/cm^. On the other hand, our resuhs are shghtly closer to those of Haritos and 
coworkers J56|, which are calculated by solving the time-independent Schrodinger equation 
by the nonperturbative many-electron, many-photon theory (MEMPT). This might be due 
to the fact that the dipole moment we use takes the core polarization effects into account. 

V. CONCLUSIONS 

We have presented an accurate and efficient method of solving the time-dependent 
Schrodinger equation. Compared to the usual FD discretization scheme of the radial co- 
ordinate, the present CWDVR method needs 3-10 fewer number of grid points and treats 
the Coulomb singularity naturally and accurately. As examples, our method has been very 
successfully applied to the multiphoton ionization dynamics of both H and H~. The total 
ionization rates estimated from the present method are in perfect agreement with other theo- 
retical results. We have also applied our method to investigate the short-time excitation and 
ionization dynamics of H in weak and strong static electric fields. The ground state survival 
probability and the ionization rate calculated using the present method are in full agree- 
ment with those results calculated by complex rotation method. Since the CWDVR treats 
the Coulomb potential accurately and needs much fewer points in r coordinate, the present 
methods opens the way to a more efficient treatment of many-electron atomic systems. 
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